<?php

/**
 * @package JAMA
 *
 *    For an m-by-n matrix A with m >= n, the singular value decomposition is
 *    an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
 *    an n-by-n orthogonal matrix V so that A = U*S*V'.
 *
 *    The singular values, sigma[$k] = S[$k][$k], are ordered so that
 *    sigma[0] >= sigma[1] >= ... >= sigma[n-1].
 *
 *    The singular value decompostion always exists, so the constructor will
 *    never fail.  The matrix condition number and the effective numerical
 *    rank can be computed from this decomposition.
 *
 * @author  Paul Meagher
 * @license PHP v3.0
 * @version 1.1
 */
class SingularValueDecomposition {

	/**
	 *    Internal storage of U.
	 *
	 * @var array
	 */
	private $U = array();

	/**
	 *    Internal storage of V.
	 *
	 * @var array
	 */
	private $V = array();

	/**
	 *    Internal storage of singular values.
	 *
	 * @var array
	 */
	private $s = array();

	/**
	 *    Row dimension.
	 *
	 * @var int
	 */
	private $m;

	/**
	 *    Column dimension.
	 *
	 * @var int
	 */
	private $n;


	/**
	 *    Construct the singular value decomposition
	 *
	 *    Derived from LINPACK code.
	 *
	 * @param $A Rectangular matrix
	 *
	 * @return Structure to access U, S and V.
	 */
	public function __construct($Arg) {

		// Initialize.
		$A = $Arg->getArrayCopy();
		$this->m = $Arg->getRowDimension();
		$this->n = $Arg->getColumnDimension();
		$nu = min($this->m, $this->n);
		$e = array();
		$work = array();
		$wantu = true;
		$wantv = true;
		$nct = min($this->m - 1, $this->n);
		$nrt = max(0, min($this->n - 2, $this->m));

		// Reduce A to bidiagonal form, storing the diagonal elements
		// in s and the super-diagonal elements in e.
		for ($k = 0; $k < max($nct, $nrt); ++$k)
		{

			if ($k < $nct)
			{
				// Compute the transformation for the k-th column and
				// place the k-th diagonal in s[$k].
				// Compute 2-norm of k-th column without under/overflow.
				$this->s[$k] = 0;
				for ($i = $k; $i < $this->m; ++$i)
				{
					$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
				}
				if ($this->s[$k] != 0.0)
				{
					if ($A[$k][$k] < 0.0)
					{
						$this->s[$k] = -$this->s[$k];
					}
					for ($i = $k; $i < $this->m; ++$i)
					{
						$A[$i][$k] /= $this->s[$k];
					}
					$A[$k][$k] += 1.0;
				}
				$this->s[$k] = -$this->s[$k];
			}

			for ($j = $k + 1; $j < $this->n; ++$j)
			{
				if (($k < $nct) & ($this->s[$k] != 0.0))
				{
					// Apply the transformation.
					$t = 0;
					for ($i = $k; $i < $this->m; ++$i)
					{
						$t += $A[$i][$k] * $A[$i][$j];
					}
					$t = -$t / $A[$k][$k];
					for ($i = $k; $i < $this->m; ++$i)
					{
						$A[$i][$j] += $t * $A[$i][$k];
					}
					// Place the k-th row of A into e for the
					// subsequent calculation of the row transformation.
					$e[$j] = $A[$k][$j];
				}
			}

			if ($wantu AND ($k < $nct))
			{
				// Place the transformation in U for subsequent back
				// multiplication.
				for ($i = $k; $i < $this->m; ++$i)
				{
					$this->U[$i][$k] = $A[$i][$k];
				}
			}

			if ($k < $nrt)
			{
				// Compute the k-th row transformation and place the
				// k-th super-diagonal in e[$k].
				// Compute 2-norm without under/overflow.
				$e[$k] = 0;
				for ($i = $k + 1; $i < $this->n; ++$i)
				{
					$e[$k] = hypo($e[$k], $e[$i]);
				}
				if ($e[$k] != 0.0)
				{
					if ($e[$k + 1] < 0.0)
					{
						$e[$k] = -$e[$k];
					}
					for ($i = $k + 1; $i < $this->n; ++$i)
					{
						$e[$i] /= $e[$k];
					}
					$e[$k + 1] += 1.0;
				}
				$e[$k] = -$e[$k];
				if (($k + 1 < $this->m) AND ($e[$k] != 0.0))
				{
					// Apply the transformation.
					for ($i = $k + 1; $i < $this->m; ++$i)
					{
						$work[$i] = 0.0;
					}
					for ($j = $k + 1; $j < $this->n; ++$j)
					{
						for ($i = $k + 1; $i < $this->m; ++$i)
						{
							$work[$i] += $e[$j] * $A[$i][$j];
						}
					}
					for ($j = $k + 1; $j < $this->n; ++$j)
					{
						$t = -$e[$j] / $e[$k + 1];
						for ($i = $k + 1; $i < $this->m; ++$i)
						{
							$A[$i][$j] += $t * $work[$i];
						}
					}
				}
				if ($wantv)
				{
					// Place the transformation in V for subsequent
					// back multiplication.
					for ($i = $k + 1; $i < $this->n; ++$i)
					{
						$this->V[$i][$k] = $e[$i];
					}
				}
			}
		}

		// Set up the final bidiagonal matrix or order p.
		$p = min($this->n, $this->m + 1);
		if ($nct < $this->n)
		{
			$this->s[$nct] = $A[$nct][$nct];
		}
		if ($this->m < $p)
		{
			$this->s[$p - 1] = 0.0;
		}
		if ($nrt + 1 < $p)
		{
			$e[$nrt] = $A[$nrt][$p - 1];
		}
		$e[$p - 1] = 0.0;
		// If required, generate U.
		if ($wantu)
		{
			for ($j = $nct; $j < $nu; ++$j)
			{
				for ($i = 0; $i < $this->m; ++$i)
				{
					$this->U[$i][$j] = 0.0;
				}
				$this->U[$j][$j] = 1.0;
			}
			for ($k = $nct - 1; $k >= 0; --$k)
			{
				if ($this->s[$k] != 0.0)
				{
					for ($j = $k + 1; $j < $nu; ++$j)
					{
						$t = 0;
						for ($i = $k; $i < $this->m; ++$i)
						{
							$t += $this->U[$i][$k] * $this->U[$i][$j];
						}
						$t = -$t / $this->U[$k][$k];
						for ($i = $k; $i < $this->m; ++$i)
						{
							$this->U[$i][$j] += $t * $this->U[$i][$k];
						}
					}
					for ($i = $k; $i < $this->m; ++$i)
					{
						$this->U[$i][$k] = -$this->U[$i][$k];
					}
					$this->U[$k][$k] = 1.0 + $this->U[$k][$k];
					for ($i = 0; $i < $k - 1; ++$i)
					{
						$this->U[$i][$k] = 0.0;
					}
				}
				else
				{
					for ($i = 0; $i < $this->m; ++$i)
					{
						$this->U[$i][$k] = 0.0;
					}
					$this->U[$k][$k] = 1.0;
				}
			}
		}

		// If required, generate V.
		if ($wantv)
		{
			for ($k = $this->n - 1; $k >= 0; --$k)
			{
				if (($k < $nrt) AND ($e[$k] != 0.0))
				{
					for ($j = $k + 1; $j < $nu; ++$j)
					{
						$t = 0;
						for ($i = $k + 1; $i < $this->n; ++$i)
						{
							$t += $this->V[$i][$k] * $this->V[$i][$j];
						}
						$t = -$t / $this->V[$k + 1][$k];
						for ($i = $k + 1; $i < $this->n; ++$i)
						{
							$this->V[$i][$j] += $t * $this->V[$i][$k];
						}
					}
				}
				for ($i = 0; $i < $this->n; ++$i)
				{
					$this->V[$i][$k] = 0.0;
				}
				$this->V[$k][$k] = 1.0;
			}
		}

		// Main iteration loop for the singular values.
		$pp = $p - 1;
		$iter = 0;
		$eps = pow(2.0, -52.0);

		while ($p > 0)
		{
			// Here is where a test for too many iterations would go.
			// This section of the program inspects for negligible
			// elements in the s and e arrays.  On completion the
			// variables kase and k are set as follows:
			// kase = 1  if s(p) and e[k-1] are negligible and k<p
			// kase = 2  if s(k) is negligible and k<p
			// kase = 3  if e[k-1] is negligible, k<p, and
			//           s(k), ..., s(p) are not negligible (qr step).
			// kase = 4  if e(p-1) is negligible (convergence).
			for ($k = $p - 2; $k >= -1; --$k)
			{
				if ($k == -1)
				{
					break;
				}
				if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k + 1])))
				{
					$e[$k] = 0.0;
					break;
				}
			}
			if ($k == $p - 2)
			{
				$kase = 4;
			}
			else
			{
				for ($ks = $p - 1; $ks >= $k; --$ks)
				{
					if ($ks == $k)
					{
						break;
					}
					$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks - 1]) : 0.);
					if (abs($this->s[$ks]) <= $eps * $t)
					{
						$this->s[$ks] = 0.0;
						break;
					}
				}
				if ($ks == $k)
				{
					$kase = 3;
				}
				else if ($ks == $p - 1)
				{
					$kase = 1;
				}
				else
				{
					$kase = 2;
					$k = $ks;
				}
			}
			++$k;

			// Perform the task indicated by kase.
			switch ($kase)
			{
				// Deflate negligible s(p).
				case 1:
					$f = $e[$p - 2];
					$e[$p - 2] = 0.0;
					for ($j = $p - 2; $j >= $k; --$j)
					{
						$t = hypo($this->s[$j], $f);
						$cs = $this->s[$j] / $t;
						$sn = $f / $t;
						$this->s[$j] = $t;
						if ($j != $k)
						{
							$f = -$sn * $e[$j - 1];
							$e[$j - 1] = $cs * $e[$j - 1];
						}
						if ($wantv)
						{
							for ($i = 0; $i < $this->n; ++$i)
							{
								$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p - 1];
								$this->V[$i][$p - 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p - 1];
								$this->V[$i][$j] = $t;
							}
						}
					}
					break;
				// Split at negligible s(k).
				case 2:
					$f = $e[$k - 1];
					$e[$k - 1] = 0.0;
					for ($j = $k; $j < $p; ++$j)
					{
						$t = hypo($this->s[$j], $f);
						$cs = $this->s[$j] / $t;
						$sn = $f / $t;
						$this->s[$j] = $t;
						$f = -$sn * $e[$j];
						$e[$j] = $cs * $e[$j];
						if ($wantu)
						{
							for ($i = 0; $i < $this->m; ++$i)
							{
								$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k - 1];
								$this->U[$i][$k - 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k - 1];
								$this->U[$i][$j] = $t;
							}
						}
					}
					break;
				// Perform one qr step.
				case 3:
					// Calculate the shift.
					$scale = max(max(max(max(abs($this->s[$p - 1]), abs($this->s[$p - 2])), abs($e[$p - 2])), abs($this->s[$k])), abs($e[$k]));
					$sp = $this->s[$p - 1] / $scale;
					$spm1 = $this->s[$p - 2] / $scale;
					$epm1 = $e[$p - 2] / $scale;
					$sk = $this->s[$k] / $scale;
					$ek = $e[$k] / $scale;
					$b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
					$c = ($sp * $epm1) * ($sp * $epm1);
					$shift = 0.0;
					if (($b != 0.0) || ($c != 0.0))
					{
						$shift = sqrt($b * $b + $c);
						if ($b < 0.0)
						{
							$shift = -$shift;
						}
						$shift = $c / ($b + $shift);
					}
					$f = ($sk + $sp) * ($sk - $sp) + $shift;
					$g = $sk * $ek;
					// Chase zeros.
					for ($j = $k; $j < $p - 1; ++$j)
					{
						$t = hypo($f, $g);
						$cs = $f / $t;
						$sn = $g / $t;
						if ($j != $k)
						{
							$e[$j - 1] = $t;
						}
						$f = $cs * $this->s[$j] + $sn * $e[$j];
						$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
						$g = $sn * $this->s[$j + 1];
						$this->s[$j + 1] = $cs * $this->s[$j + 1];
						if ($wantv)
						{
							for ($i = 0; $i < $this->n; ++$i)
							{
								$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1];
								$this->V[$i][$j + 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j + 1];
								$this->V[$i][$j] = $t;
							}
						}
						$t = hypo($f, $g);
						$cs = $f / $t;
						$sn = $g / $t;
						$this->s[$j] = $t;
						$f = $cs * $e[$j] + $sn * $this->s[$j + 1];
						$this->s[$j + 1] = -$sn * $e[$j] + $cs * $this->s[$j + 1];
						$g = $sn * $e[$j + 1];
						$e[$j + 1] = $cs * $e[$j + 1];
						if ($wantu && ($j < $this->m - 1))
						{
							for ($i = 0; $i < $this->m; ++$i)
							{
								$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j + 1];
								$this->U[$i][$j + 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j + 1];
								$this->U[$i][$j] = $t;
							}
						}
					}
					$e[$p - 2] = $f;
					$iter = $iter + 1;
					break;
				// Convergence.
				case 4:
					// Make the singular values positive.
					if ($this->s[$k] <= 0.0)
					{
						$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
						if ($wantv)
						{
							for ($i = 0; $i <= $pp; ++$i)
							{
								$this->V[$i][$k] = -$this->V[$i][$k];
							}
						}
					}
					// Order the singular values.
					while ($k < $pp)
					{
						if ($this->s[$k] >= $this->s[$k + 1])
						{
							break;
						}
						$t = $this->s[$k];
						$this->s[$k] = $this->s[$k + 1];
						$this->s[$k + 1] = $t;
						if ($wantv AND ($k < $this->n - 1))
						{
							for ($i = 0; $i < $this->n; ++$i)
							{
								$t = $this->V[$i][$k + 1];
								$this->V[$i][$k + 1] = $this->V[$i][$k];
								$this->V[$i][$k] = $t;
							}
						}
						if ($wantu AND ($k < $this->m - 1))
						{
							for ($i = 0; $i < $this->m; ++$i)
							{
								$t = $this->U[$i][$k + 1];
								$this->U[$i][$k + 1] = $this->U[$i][$k];
								$this->U[$i][$k] = $t;
							}
						}
						++$k;
					}
					$iter = 0;
					--$p;
					break;
			} // end switch
		} // end while

	} // end constructor


	/**
	 *    Return the left singular vectors
	 *
	 * @access public
	 * @return U
	 */
	public function getU() {
		return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
	}


	/**
	 *    Return the right singular vectors
	 *
	 * @access public
	 * @return V
	 */
	public function getV() {
		return new Matrix($this->V);
	}


	/**
	 *    Return the one-dimensional array of singular values
	 *
	 * @access public
	 * @return diagonal of S.
	 */
	public function getSingularValues() {
		return $this->s;
	}


	/**
	 *    Return the diagonal matrix of singular values
	 *
	 * @access public
	 * @return S
	 */
	public function getS() {
		for ($i = 0; $i < $this->n; ++$i)
		{
			for ($j = 0; $j < $this->n; ++$j)
			{
				$S[$i][$j] = 0.0;
			}
			$S[$i][$i] = $this->s[$i];
		}

		return new Matrix($S);
	}


	/**
	 *    Two norm
	 *
	 * @access public
	 * @return max(S)
	 */
	public function norm2() {
		return $this->s[0];
	}


	/**
	 *    Two norm condition number
	 *
	 * @access public
	 * @return max(S)/min(S)
	 */
	public function cond() {
		return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
	}


	/**
	 *    Effective numerical matrix rank
	 *
	 * @access public
	 * @return Number of nonnegligible singular values.
	 */
	public function rank() {
		$eps = pow(2.0, -52.0);
		$tol = max($this->m, $this->n) * $this->s[0] * $eps;
		$r = 0;
		for ($i = 0; $i < count($this->s); ++$i)
		{
			if ($this->s[$i] > $tol)
			{
				++$r;
			}
		}

		return $r;
	}

} //	class SingularValueDecomposition
